%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                         %
%                             Main Program                                %
%                                                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% ---------------------------- Description -------------------------------- 

    % This program calibrate the model to the initial balanced growth path
    % and the ending balanced growth path. The decomposition exercise
    % identifies the contributions of each potential channels. The program
    % generates Table 2-9 in the paper.

clc
clear all

global eta lambda zeta
global chi rho gamma mu iota
global phi nu theta
global epsilon beta delta delta_p
global Q alphah alphal mh ml zbar
global ss

options=optimoptions('fsolve','Display','off');
options2=optimoptions('fmincon','OptimalityTolerance',1E-8,'Display','off');


%--------------------------------------------------------------------------
%   Parameters Values from a Priori Information or Direct Estimation      %
%--------------------------------------------------------------------------

eta=0.85*(1/3);        % Capital Share
lambda=0.85*(2/3);     % Labor Share
zeta=1-eta-lambda;     % Profit Share

epsilon=2;             % CRRA Parameter
delta=0.069;           % Depreciation Rate
delta_p=0.94;          % Patent Survival Rate

r=1/(1+0.06);
g=(1+0.0213)^((zeta+lambda)/zeta);
beta=r*g^(epsilon*zeta/(zeta+lambda)); % Discount Factor

Q=[0.872, 1-0.872; 0.017, 1-0.017];    % Type Transition Matrix
ptt=(1/2)*ones(1,2)';
tol=10^-6;

% Find the stationary distribution of idiosyncratic shocks
for j=2:1000
    pptt=Q'*ptt;
    if max(abs(pptt-ptt))<tol
        ptt_star=ptt;
        break
    end
    ptt=pptt;
end

alphah=ptt_star(1);    % Share of High Type
alphal=ptt_star(2);    % Share of Low Type
ml=exp(-0.355);        % Prod. Ability of High Type
mh=exp(3.196);         % Prod. Ability of Low Type

chi=1;                 % R&D Cost, Scale, Normalization
zbar=1;                % Average innovation level, Normalization
nu=0.7;                % Matching Function, Exponent

%--------------------------------------------------------------------------
%                        Setting Data Targets                             %
%--------------------------------------------------------------------------

% Initial BGP
g_target=(1+0.0213)^((zeta+lambda)/zeta);
rd_sales_h=0.0362; 
rd_sales_l=0.0283;
range_h=11.81;     
range_l=1.924;
range_avg=alphah*range_h+alphal*range_l;
trans_target=0.2320;

% Ending BGP
g_target2=(1+0.0222)^((zeta+lambda)/zeta);
rd_sales_h2=0.0315;   
rd_sales_l2=0.0671;   
range_h2=6.307;        
range_l2=1.609;
range_avg2=alphah*range_h2+alphal*range_l2;
trans_target2=0.3704;


%--------------------------------------------------------------------------
%                     Calibrating to the Initial BGP                      %
%--------------------------------------------------------------------------

ss=0.05; % R&D Tax Credit in 1981-1985

% Calibrate key parameters

par_vec=@(vec)par_solver_policy(vec,g_target,rd_sales_h,rd_sales_l,range_h,range_l,trans_target);
vec_guess=[2 ,   0.0002  , 1.5  ,  0.8  ,  0.05  ,  0.2]; 

lb=[0,0,0,0,0,0];
ub=[1E6,1E6,1E6,1E6,1E6,1];
[par_opt,V_opt,exitflag]=fmincon(par_vec,vec_guess,[],[],[],[],[],[],[],options2);

if exitflag <=0
    disp('Parameters Not Solved')
end

iota=par_opt(1);  % Management Cost, Elasticity
mu=par_opt(2);    % Management Cost, Scale
gamma=par_opt(3); % Step Size of Innovation
rho=par_opt(4);   % 1+rho: R&D Cost Elasticity
phi=par_opt(5);   % Matching Function, Scale
theta=par_opt(6); % Bargaining Power of the Patent Seller

% Compute moments in the initial BGP
guess=[0.9,0.1,11,2,1,6,1,2];
geq_vec=@(input)geq_solver_policy(input,iota,mu,gamma,rho,phi,theta);
[answer,value,exitflag] = fsolve(geq_vec, guess, options);
if exitflag ~= 1
    disp('First-Order Optimality is not Achieved')
end

ih=answer(1);
il=answer(2);
omegah=answer(3);
omegal=answer(4);
a=answer(5);
v1h=answer(6);
v1l=answer(7);
ns=answer(8);

nbh=alphah*(1-ih*X(omegah));
nbl=alphal*(1-il*X(omegal));
nb=nbh+nbl;
sigmah=nbh/nb;
sigmal=nbl/nb;
mb=phi*(ns/nb)^nu;
ms=phi*(nb/ns)^(1-nu);

eh=(alphah*mh*a+alphal*ml)*(alphah*Q(1,1)+alphal*Q(2,1))/((alphah*mh+alphal*ml)*(alphah*Q(1,1)*a+alphal*Q(2,1)));
el=(alphah*mh*a+alphal*ml)*(alphal*Q(2,2)+alphah*Q(1,2))/((alphah*mh+alphal*ml)*(alphal*Q(2,2)+alphah*Q(1,2)*a));
gh=1+gamma*eh*(ih*X(omegah)+(1-ih*X(omegah))*mb);
gl=1+gamma*el*(il*X(omegal)+(1-il*X(omegal))*mb);

g_avg=gh*(alphah*Q(1,1)+alphal*Q(2,1)/a)/(alphah*Q(1,1)+alphal*Q(2,1));
rm=beta/(g_avg^(epsilon*zeta/(lambda+zeta)));
rt=1/rm-1+delta;
w=lambda*((eta/rt)^(eta/(zeta+lambda)))*((alphah*mh+alphal*ml)*g_avg*zbar)^(zeta/(lambda+zeta));

o=(ms*theta*(sigmah*v1h+sigmal*v1l)*gamma)/(1-(1-ms*theta)*rm*delta_p*g_avg^(zeta/(lambda+zeta)));

zlbar=(alphah*mh+alphal*ml)*zbar/(alphah*mh*a+alphal*ml);
zhbar=a*zlbar;

buy_share_h=(1-ih*X(omegah))*mb/ih;                
buy_share_l=(1-il*X(omegal))*mb/il;

sell_share_h=ih*(1-X(omegah))*ms/ih;
sell_share_l=il*(1-X(omegal))*ms/il;
    
rd_exp_h=(1-ss)*chi*(zbar^(zeta/(zeta+lambda)))*(ih^(1+rho))/(1+rho);
rd_exp_l=(1-ss)*(zbar^(zeta/(zeta+lambda)))*(il^(1+rho))/(1+rho);

zhr=(alphah*Q(1,1)*zhbar+alphal*Q(2,1)*zlbar)/(alphah*Q(1,1)+alphal*Q(2,1));
zlr=(alphah*Q(1,2)*zhbar+alphal*Q(2,2)*zlbar)/(alphah*Q(1,2)+alphal*Q(2,2));

Yh=mh*gh*zhr*((eta/rt)^(eta/zeta))*((lambda/w)^(lambda/zeta));
Yl=ml*gl*zlr*((eta/rt)^(eta/zeta))*((lambda/w)^(lambda/zeta));

rd_Y_h=rd_exp_h/Yh;
rd_Y_l=rd_exp_l/Yl;

within_scope_h_initial=X(omegah);
within_scope_l_initial=X(omegal);

patent_count=alphah*ih+alphal*il;
trans_share=(alphah*ih*(1-X(omegah))+alphal*il*(1-X(omegal)))*(ms*(1-((1-ms)*delta_p)^10))/((1-(1-ms)*delta_p)*patent_count);

if v1h*gamma >= o && v1l*gamma>=o
    disp('Both types keep patents');
elseif v1h*gamma > o && v1l*gamma<o
    disp('High type keeps patents; low type sells patents');
else
    disp('Both types sell patents');
end


%--------------------------------------------------------------------------
%                            Generate Tables                              %
%--------------------------------------------------------------------------

phi_initial=phi; theta_initial=theta; mu_initial=mu; iota_initial=iota; gamma_initial=gamma; rho_initial=rho; ss_initial=ss;
omegah_initial=omegah; omegal_initial=omegal; omega_avg_initial=alphah*omegah+alphal*omegal; rd_Y_h_initial=rd_Y_h; rd_Y_l_initial=rd_Y_l; trans_share_initial=trans_share; g_avg_initial=g_avg;

% Create Table 2: Parameter Values from a Priori Information
Parameter = {
    'η';'λ';'ε';'β';'δ';'α_H';'α_L';'χ';'σ'};

Description = {
    'Capital Share';'Labor Share';'CRRA Parameter';'Discount Factor';'Depreciation Rate';'Share of High Type';'Share of Low Type';'R&D Cost Scale';'R&D Tax Credit Rate'};

Value = [eta; lambda; epsilon; beta; delta; alphah; alphal; chi; ss];

Identification = {
    'Guner et al. (2008)';
    'Guner et al. (2008)';
    'Standard';
    'Interest Rate Target';
    'NIPA';
    'Imposed';
    'Imposed';
    'Normalization';
    'Akcigit et al. (2018)'
};

Table2 = table(Parameter, Description, Value, Identification, ...
    'VariableNames', {'Parameter', 'Description', 'Value', 'Identification'});
disp(Table2);
writetable(Table2, "C:\Users\Yueyuan Ma\Box\Specialization\Replication_Files\Tables\Table2.xlsx");


% Create Table 3: Parameter Values from Direct Estimations
Parameter = {
    'mH';'mL';'Qmm';'nu';'X(omega)'};

Description = {
    'Prod. Ability of High Type';
    'Prod. Ability of Low Type';
    'Type Transition Matrix';
    'Matching Function, Exponent';
    'Within-scope Probability'
};

Value = { mh; ml;  Q;  nu;  "exp(-4.443)*|omega|^0.7643;"};

Identification = {
    'Regression';
    'Regression';
    'MLE';
    'Regression';
    'Regression'
};

Table3 = table(Parameter, Description, Value, Identification, ...
    'VariableNames', {'Parameter', 'Description', 'Value', 'Identification'});
disp(Table3);
writetable(Table3, "C:\Users\Yueyuan Ma\Box\Specialization\Replication_Files\Tables\Table3.xlsx");


% Create Table 4: Parameter Values from the Minimum Distance Estimation
Parameter = {'γ'; '1+ρ'; 'θ'; 'µ'; '1+ι'; 'ϕ'};

Description = {
    'Step Size of Innovation';
    'R&D Cost Elasticity';
    'Bargaining Power';
    'Management Cost, Scale';
    'Management Cost, Elasticity';
    'Matching Function, Scale'
};

Value = [gamma; 1+rho; theta; mu; 1+iota; phi];

Identification = {
    'Growth Rate';
    'R&D Cost/Sales';
    'Ratio (H and L)';
    'Avg. Number of Industries (H and L)';
    'Avg. Number of Industries (H and L)';
    'Patent Traded Share'
};


Table4 = table(Parameter, Description, Value, Identification, ...
    'VariableNames', {'Parameter', 'Description', 'Value', 'Identification'});
disp(Table4);
writetable(Table4, "C:\Users\Yueyuan Ma\Box\Specialization\Replication_Files\Tables\Table4.xlsx");

% Create Table 5: Model Fit for Key Moments in the Initial Balanced Growth Path

Targets = {
    'Economic Growth Rate (1981-1985)';
    'R&D Cost/Sales of H Firms (1981-1985)';
    'R&D Cost/Sales of L Firms (1981-1985)';
    'Avg. Number of Industries of H Firms (1981-1985)';
    'Avg. Number of Industries of L Firms (1981-1985)';
    'The Share of Patents Transacted (1981-1985)'
};

Data = {g_target^(zeta/(zeta+lambda))-1; rd_sales_h; rd_sales_l; range_h; range_l; trans_target};

Model = {g_avg^(zeta/(zeta+lambda))-1; rd_Y_h; rd_Y_l; omegah; omegal; trans_share};

Table5 = table(Targets, Data, Model, ...
    'VariableNames', {'Targets', 'Data', 'Model'});
disp(Table5);
writetable(Table5, "C:\Users\Yueyuan Ma\Box\Specialization\Replication_Files\Tables\Table5.xlsx");


%--------------------------------------------------------------------------
%                     Calibrating to the Ending BGP                      %
%--------------------------------------------------------------------------

ss=0.24; % R&D Tax Credit in 1996-2000

% Calibrate key parameters

par_vec=@(vec)par_solver_policy2(vec,g_target2,rd_sales_h2,rd_sales_l2,range_h2,range_l2,trans_target2);
vec_guess=[2 ,   0.0002  , 1.5  ,  0.8  ,  0.1  ,  0.2]; 

lb=[0,0,0,0,0,0];
ub=[1E6,1E6,1E6,1E6,1E6,1];
[par_opt,V_opt,exitflag]=fmincon(par_vec,vec_guess,[],[],[],[],[],[],[],options2);

if exitflag <=0
    disp('Parameters Not Solved')
end

iota=par_opt(1);
mu=par_opt(2);
gamma=par_opt(3);
rho=par_opt(4);
phi=par_opt(5);
theta=par_opt(6);

% Compute Moments in the ending BGP
guess=[0.7,0.1,6,1.5,1,5,1,0.4];
geq_vec=@(input)geq_solver_policy(input,iota,mu,gamma,rho,phi,theta);
[answer,value,exitflag] = fsolve(geq_vec, guess, options);
if exitflag ~= 1
    disp('First-Order Optimality is not Achieved')
end

ih=answer(1);
il=answer(2);
omegah=answer(3);
omegal=answer(4);
a=answer(5);
v1h=answer(6);
v1l=answer(7);
ns=answer(8);

nbh=alphah*(1-ih*X(omegah));
nbl=alphal*(1-il*X(omegal));
nb=nbh+nbl;
sigmah=nbh/nb;
sigmal=nbl/nb;
mb=phi*(ns/nb)^nu;
ms=phi*(nb/ns)^(1-nu);

eh=(alphah*mh*a+alphal*ml)*(alphah*Q(1,1)+alphal*Q(2,1))/((alphah*mh+alphal*ml)*(alphah*Q(1,1)*a+alphal*Q(2,1)));
el=(alphah*mh*a+alphal*ml)*(alphal*Q(2,2)+alphah*Q(1,2))/((alphah*mh+alphal*ml)*(alphal*Q(2,2)+alphah*Q(1,2)*a));
gh=1+gamma*eh*(ih*X(omegah)+(1-ih*X(omegah))*mb);
gl=1+gamma*el*(il*X(omegal)+(1-il*X(omegal))*mb);

g_avg=gh*(alphah*Q(1,1)+alphal*Q(2,1)/a)/(alphah*Q(1,1)+alphal*Q(2,1));
rm=beta/(g_avg^(epsilon*zeta/(lambda+zeta)));
rt=1/rm-1+delta;
w=lambda*((eta/rt)^(eta/(zeta+lambda)))*((alphah*mh+alphal*ml)*g_avg*zbar)^(zeta/(lambda+zeta));

o=(ms*theta*(sigmah*v1h+sigmal*v1l)*gamma)/(1-(1-ms*theta)*rm*delta_p*g_avg^(zeta/(lambda+zeta)));

zlbar=(alphah*mh+alphal*ml)*zbar/(alphah*mh*a+alphal*ml);
zhbar=a*zlbar;

buy_share_h=(1-ih*X(omegah))*mb/ih;               
buy_share_l=(1-il*X(omegal))*mb/il;

sell_share_h=ih*(1-X(omegah))*ms/ih;
sell_share_l=il*(1-X(omegal))*ms/il;
    
rd_exp_h=(1-ss)*chi*(zbar^(zeta/(zeta+lambda)))*(ih^(1+rho))/(1+rho);
rd_exp_l=(1-ss)*chi*(zbar^(zeta/(zeta+lambda)))*(il^(1+rho))/(1+rho);

zhr=(alphah*Q(1,1)*zhbar+alphal*Q(2,1)*zlbar)/(alphah*Q(1,1)+alphal*Q(2,1));
zlr=(alphah*Q(1,2)*zhbar+alphal*Q(2,2)*zlbar)/(alphah*Q(1,2)+alphal*Q(2,2));

Yh=mh*gh*zhr*((eta/rt)^(eta/zeta))*((lambda/w)^(lambda/zeta));
Yl=ml*gl*zlr*((eta/rt)^(eta/zeta))*((lambda/w)^(lambda/zeta));

rd_Y_h=rd_exp_h/Yh;
rd_Y_l=rd_exp_l/Yl;

within_scope_h_ending=X(omegah);
within_scope_l_ending=X(omegal);

patent_count=alphah*ih+alphal*il;
trans_share=(alphah*ih*(1-X(omegah))+alphal*il*(1-X(omegal)))*(ms*(1-((1-ms)*delta_p)^10))/((1-(1-ms)*delta_p)*patent_count);

if v1h*gamma >= o && v1l*gamma>=o
    disp('Both types keep patents');
elseif v1h*gamma > o && v1l*gamma<o
    disp('High type keeps patents; low type sells patents');
else
    disp('Both types sell patents');
end

%--------------------------------------------------------------------------
%                            Generate Tables                              %
%--------------------------------------------------------------------------

phi_ending=phi; theta_ending=theta; mu_ending=mu; iota_ending=iota; gamma_ending=gamma; rho_ending=rho; ss_ending=ss;

% Create Table 6: Model Fit for Key Moments in the Ending Balanced Growth Path

Targets = {
    'Economic Growth Rate (1996-2000)';
    'R&D Cost/Sales of H Firms (1996-2000)';
    'R&D Cost/Sales of L Firms (1996-2000)';
    'Avg. Number of Industries of H Firms (1996-2000)';
    'Avg. Number of Industries of L Firms (1996-2000)';
    'The Share of Patents Transacted (2000)'
};

Data = {g_target2^(zeta/(zeta+lambda))-1; rd_sales_h2; rd_sales_l2; range_h2; range_l2; trans_target2};

Model = {g_avg^(zeta/(zeta+lambda))-1; rd_Y_h; rd_Y_l; omegah; omegal; trans_share};

Table6 = table(Targets, Data, Model, ...
    'VariableNames', {'Targets', 'Data', 'Model'});
disp(Table6);
writetable(Table6, "C:\Users\Yueyuan Ma\Box\Specialization\Replication_Files\Tables\Table6.xlsx");


% Create Table 7: Model Fit for Untargeted Moments

Moments = {
    'Within-scope Prob. of H Firms (1981-1985)';
    'Within-scope Prob. of L Firms (1981-1985)';
    'Within-scope Prob. of H Firms (1996-2000)';
    'Within-scope Prob. of L Firms (1996-2000)'
};

Data = {
    '6.65%';   % H Firms 1981–85
    '2.92%';   % L Firms 1981–85
    '3.79%';   % H Firms 1996–00
    '2.25%'    % L Firms 1996–00
};

Model = {
    within_scope_h_initial;   % H Firms 1981–85
    within_scope_l_initial;   % L Firms 1981–85
    within_scope_h_ending;    % H Firms 1996–00
    within_scope_l_ending     % L Firms 1996–00
};

Table7 = table(Moments, Data, Model, ...
    'VariableNames', {'Moments', 'Data', 'Model'});
disp(Table7);
writetable(Table7, "C:\Users\Yueyuan Ma\Box\Specialization\Replication_Files\Tables\Table7.xlsx");


% Create Table 8: Changes of Parameter Values
Parameter = {'ϕ'; 'θ'; 'µ'; '1 + ι'; 'γ'; '1 + ρ'};

Old_BGP = [phi_initial; theta_initial; mu_initial; 1+iota_initial; gamma_initial; 1+rho_initial];

New_BGP = [phi_ending; theta_ending; mu_ending; 1+iota_ending; gamma_ending; 1+rho_ending];

Interpretation = {
    'Matching efficiency increase';
    'Sellers’’ bargaining power increase';
    'Higher costs of expanding scope';
    'More decreasing return to scope';
    'Fall in R&D efficiency';
    'Fall in R&D efficiency'
};


Table8 = table(Parameter, Old_BGP, New_BGP, Interpretation, ...
    'VariableNames', {'Parameter', 'Old_BGP', 'New_BGP', 'Interpretation'});
disp(Table8);
writetable(Table8, "C:\Users\Yueyuan Ma\Box\Specialization\Replication_Files\Tables\Table8.xlsx");

% Create Table 9: Effects of Key Parameters

Decomposition
